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Abstract 

This  paper  describes  a  3D  thermal  modeling  by  a  nodes  network  model  for  two  PEMFC  of  150  and  500  W  (respectively,  3  and  20  cells).  Modeling 
are  realized  for  each  case  for  one  cell  before  to  be  integrated  on  all  the  stack.  Absolute  temperatures  of  H2,  air  and  water  channels  are  used  as 
Dirichlet  conditions.  Temperatures  of  external  surfaces  are  obtained  thanks  to  an  infrared  thermographic  camera.  Final  external  heat  fluxes  are 
deduced  from  the  integrated  model. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Operating  conditions  of  a  fuel  cell  widely  depend  on  the 
thermal  management.  It  is  used  to  control  the  cooling  system,  to 
maintain  a  good  hygrometry  level  in  the  fuel  cell  and  to  optimize 
the  global  efficiency  of  the  system.  Some  studies  concerning 
the  solid  structures  have  been  made  by  considering  isothermal 
cells.  Then,  the  computations  of  the  temperatures  are  obtained 
in  the  fluids  of  the  ducts  [1].  Other  studies  are  made  where  the 
distribution  of  temperatures  is  computed  on  the  working  plates 
[2],  but  the  gradients  of  temperature  through  the  layers  (MEA) 
are  not  taken  into  account,  then  the  heat  transfer  can  only  be 
estimated  along  the  channels.  These  studies  are  realized  with 
two  kinds  of  boundary  conditions:  free  convection  or  with  a 
water  circulation  on  the  external  faces  (forced  convection). 
3D  thermal  computations  are  done  considering  the  imposed 
temperatures  of  the  working  plates  as  Dirichlet  boundary 
conditions  [3,4].  The  main  object  of  this  work  is  to  study  the 
influence  of  temperature  on  physical  parameters  as  for  example, 
the  hygrometry  rate  of  the  membrane.  The  cooling  system  is  not 
taken  into  account.  All  these  works  focused  on  different  parts  of 
one  cell.  This  latter  is  not  considered  as  an  entire  unit.  Different 
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ways  where  the  cell  evacuates  the  produced  heat  power  are 
studied  to  establish  a  thermal  management.  So  the  thermal 
flux  paths  must  be  known  and  particularly  the  part  of  the  heat 
which  must  be  removed  by  the  cooling  system  or  transferred 
by  conduction-convection  across  the  faces  of  the  stack.  It  is 
necessary  to  know  the  ability  of  the  heat  fluxes  to  cross  through 
the  cells  and  to  quantify  the  part  concentrated  in  the  bipolar 
plates  and  dissipated  on  the  external  surfaces.  This  study  uses  a 
3D  thermal  nodes  network  to  compute  the  temperatures  and  the 
flux  distribution  in  the  cell.  Various  boundary  conditions  are 
applied,  cyclic  or  imposed  flux  to  study  the  thermal  operating 
mode  between  two  adjoining  cells  or  the  heat  transfer  along  the 
length  of  the  stack.  The  distribution  of  temperatures  is  obtained 
for  different  current  densities  for  two  stacks  respectively  of 
150  W  (3  cells)  and  500  W  (20  cells).  The  surface  temperature  is 
supposed  as  constant  and  defined  by  the  cooling  system  whose 
temperature  is  controlled  at  the  output  and  the  heat  transfer 
through  the  faces  is  also  constant.  The  heat  flux  extracted  by 
the  cooling  system  is  therefore  greatly  depending  on  the  current 
density.  Comparisons  with  experiments  are  presented. 

2.  Thermal  modeling  by  3D  nodal  network 

The  modeling  of  the  solid  structure  of  one  cell 
(14  cm  x  14  cm  x  0.3  cm)  is  carried  out  with  a  nodal  network 
consisted  of  172  nodes  (temperatures  of  the  wall  are  measured 
and  used  as  boundary  condition)  or  236  nodes  (temperatures  of 
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Fig.  1.  Elements  of  the  nodal  network. 


the  wall  are  computed)  of  volume  and  surface.  The  nodes  rep¬ 
resent  heat  and  mass  transfer  between  fluid  and  solid  interfaces 
[5].  Fig.  1  shows  a  typical  node  with  its  thermal  conductances, 
each  of  them  can  represent  radiation,  convection,  conduction 
and/or  mass  transfer  with  surrounding  nodes. 

In  the  modeling,  radiation  is  not  taken  into  account  con¬ 
sidering  the  low  temperature  levels.  Heat  power  sources  are 
electrochemical  reactions  of  the  production  of  H2O  in  the  cath¬ 
ode  channel.  Eq.  (1)  is  the  thermal  balance  for  the  node 
given  for  1-172  (236  when  including  nodes  of  the  fluids).  In 
this  equation,  Gij  represents  a  “conductance”  in  WK-1  which 
is  the  inverse  of  the  thermal  resistance: 

c'^7  =  ZG',/(F  -  ^  +  Q>  (D 

j 


Fig.  2.  Description  of  the  cooling  system. 


Table  1 

Hydraulic  diameters  for  the  cooling  system 


Main  duct  (mm) 

7.1 

Secondary  duct  (mm) 

3.5 

Cell  duct  (mm) 

0.8 

With  the  different  expressions  of  conductances  given  by  Eqs. 
(2a)-(2c): 


Si  i 

— 

T 

^hl 

Gij 

(2a) 

Gij 

~  hi,jSi,j 

(2b) 

Gij 

=  fncp 

(2c) 

is  the  thermal  conductivity  of  the  material  Lq  and  Sq  are 
respectively  the  distance  and  the  heat  exchange  surface  between 
nodes  i  and  j.  Hq  the  coefficient  of  convection  at  the  interface  Sjj, 
m  the  mass  flow  (kg  s-1)  and  cp  the  specific  heat  (Jkg-1  K-1). 

The  final  system  to  solve  also  consists  of  172  or  236  analytical 
equations  [6],  respectively  for  the  3  and  20  cells.  The  system  is 
then  transformed  into  an  equation  composed  of  matrix  and  vec¬ 
tors  and  solved  with  a  Gear  method  thanks  to  the  mathematical 
software  MATLAB. 

The  nodal  modeling  is  carried  out  for  only  one  cell  of  the 
PEMFC.  Internal  Dirichlet  boundary  conditions  are  obtained 
thanks  to  the  input  temperatures  of  the  water  cooling  system  ( K 
thermocouple  at  the  entrance  of  the  stack). 

Fig.  2  represents  the  simplified  cooling  system  which  uses 
demineralized  water.  As  shown  in  this  figure,  modeling  of  the 
channels  is  significantly  simplified  with  only  three  ducts  instead 
of  the  real  geometry.  This  choice  is  justified  by  the  compromise 
between  the  satisfying  accuracy  on  the  temperature  distribu¬ 
tion  and  the  shorter  duration  of  computation  determined  by  the 
number  of  nodes.  However,  the  really  mass  fluxes  are  taken  into 
account  and  adapted  transfer  surfaces  are  considered  with  equiv¬ 
alent  hydraulic  diameters  given  in  Table  1 . 

3.  Model  implementation 

Fig.  3  describes  the  developed  model.  The  modeled  water 
channels  are  represented  as  channels  1-3.  These  nodal  “con- 
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Fig.  3.  Different  parts  of  the  modeling  of  the  cooling  system. 


ductances”  characterize  heat  transfer  by  convection  between  the 
fluid  and  the  surrounding  solid  surfaces  so  that  the  mass  transfer 
exists  along  the  duct.  Conductance  4  is  the  Dirichlet  condition 
imposed  at  the  input  of  the  cooling  channel  at  the  entrance  of 
the  stack.  The  input  temperature  of  the  water  is  represented  with 
the  conductance  12. 

Its  value  is  calculated  and  compared  to  the  experimental  value 
obtained  thanks  to  a  thermocouple  of  the  test  bench.  Conver¬ 
gence  between  the  simulated  and  the  experimental  temperature 
is  obtained  by  modifications  of  the  leak  conductances  1 1  and/or 
5.  It  represents  the  external  heat  losses  of  the  modeled  cell.  Con¬ 
ductances  7  and  10  represent  the  heat  fluxes  of  the  studied  cell 
with  the  surrounding  air.  Finally,  conductances  8  and  10  repre¬ 
sent  the  thermal  influence  of  the  neighboring  cells.  They  permit 
to  consider  the  thermal  sources  and/or  the  conduction  losses  in 
the  third  dimension  (axis  O ,  z). 


4.  Boundary  conditions 

Three  configurations  of  modeling  have  been  considered  for 
the  studied  cell.  Differences  come  from  different  applications  of 
boundary  conditions.  First,  applied  boundary  conditions  have 
to  represent  the  heat  fluxes  (B.C.  First  kind)  between  different 
associated  cells.  The  other  applied  boundary  conditions  on  the 
external  surfaces  are  first  kind  boundary  conditions  of  Dirichlet 
(absolute  temperature)  [7] .  The  knowledge  (or  not)  of  these  lasts 
is  at  the  origin  of  a  two  different  configurations. 


•  Configuration  1 .  Thermal  flux  lost  on  one  interface  is  gained 
on  the  opposite  one.  The  absolute  temperatures  on  the  external 
surface  are  known. 

•  Configuration  2.  Thermal  flux  lost  on  one  interface  is  gained 
on  the  opposite  one  but  the  temperatures  on  the  external  sur¬ 
faces  are  unknown  and  finally  obtained  thanks  to  an  iterative 
process. 

•  Configuration  3.  Thermal  flux  which  is  transferred  on  the 
first  interface  can  be  different  from  the  one  transferred  on  the 
opposite  interface.  Both  values  are  algebraic  expressions.  The 
external  temperatures  are  obtained  as  for  configuration  2. 


Absolute  temperatures  (Dirichlet  conditions)  on  the  external 
surfaces  are  obtained  thanks  to  an  infrared  short  waves  cam¬ 
era  (case  1  of  chapter  results).  The  emissivities  of  the  external 
surfaces  are  obtained  by  comparison  and  convergence  of  the 
temperature  obtained  with  a  thermocouple  and  the  other  one 
obtained  by  radiation.  The  convergence  is  obtained  for  an  emis- 
sivity  of  all  surfaces  at  0.94.  These  temperatures  can  be  used 
as  boundary  conditions.  It  is  possible,  as  done  in  this  study,  to 
used  them  only  to  valid  the  model  by  comparison  of  the  final 
surface  computed  temperatures  with  the  temperatures  obtained 
by  radiation  measurement. 

5.  Equations 

To  define  the  conductances  and  the  thermal  resistances,  it 
is  necessary  to  know  all  internal  power  sources.  In  that  way, 
chemical  energy  released  by  consummation  of  hydrogen  is  given 
by  Eq.  (3): 

/ 

^chemical  —  AHy[2Og——  (3) 

Zr 

Eq.  (4)  gives  the  electric  energy  supplied  to  the  load: 


^electrical  —  UI  (4) 

Finally,  internal  thermal  power  source  is  equal  to  the  differ¬ 
ence  between  (3)  and  (4)  is  shown  by  Eq.  (5): 

^thermal  —  ^chemical  ^electrical  (5) 

Thermal  conductivities  used  for  the  modeling  are  given  in 
Table  2. 

Different  correlations  are  used  to  define  the  internal  thermal 
convection  coefficients.  The  flow  in  the  ducts  is  supposed  fully 
developed.  Different  geometries  are  used  by  the  manufacturer  of 
the  bipolar.  The  geometries  in  the  main  and  secondary  cooling 
ducts  are  circular  and  rectangular  for  the  ducts  in  the  cell.  In  the 
two  first  configurations,  the  Colburn  correlation  [8]  is  used.  It  is 
expressed  thanks  to  Eq.  (6)  where  the  flow  is  turbulent: 

Nu  =  0.023  Re0  &  Prl/ 3  (6) 

where  the  Reynolds  number  is  defined  as 

UmX 

Re=—  (7) 

v 

and  the  Prandtl  number  expressed  as: 

Pr=l  (8) 

Pr  is  calculated  at  the  arithmetic  average  temperature  of  water 

between  the  inlet  and  the  outlet  of  the  cooling  system.  The  ducts 
in  channels  1-3  are  rectangular.  Mass  flow  is  considered  as 
laminar  flow  (small  value  of  Re).  In  that  case,  correlation  (9) 


Table  2 

Thermal  conductivities  of  the  materials  of  the  cell 

Bipolar  plates  (carbon),  Xc  (W  m_  1  K- 1 )  20 

Diffusion  layer,  Zdiff  (W  m_  1  K- 1 )  0.1 

Membrane,  Zmem  (W  m“ 1  K-1)  0.6 
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permitting  to  determine  the  Nusselt  number  issued  from  Shah 
and  London  [9]  is  used: 

Nu  =  7.541(1  -  2.610£  +  4.970£2  -  5.1 19/  +  2.702£4 

—  0.548£5)  (9) 

with  a  shape  coefficient  f  of  0.7  defined  as  the  ratio  of  the  length 
of  the  cell  along  the  axes  O ,  x  and  O ,  y  (Fig.  2). 

Nusselt  numbers  permit  to  calculate  the  heat  transfer  coef¬ 
ficients  by  free  convection  on  the  external  surface  of  the  cell. 
Nu  depending  on  the  surface  temperatures,  an  iterative  method 
has  then  been  used  in  configurations  2  and  3.  At  each  step, 
it  has  been  calculated  with  the  new  value  of  the  temperature 
issued  from  the  precedent  step  of  the  computation.  The  conver¬ 
gence  is  accepted  for  a  0.1%  variation  of  temperature.  Then,  for 
4.2  x  106  <Ra<l .6  x  106  [10],  Nu  is  given  by  Eq.  (10)  for  the 
horizontal  surface  at  the  top  of  the  cell: 

Nu  =  0.27  Ral/A  (10) 


with 


Ra  =  Gr  Pr  = 


<g/3|Tp  T0 o\x^ 

av 


(11) 


rp  is  the  wall  temperature  and  T0 0  is  the  temperature  of  the 
external  surrounding: 
p  is  defined  as: 


fd  +  tc 


oo 


(12) 


For  the  horizontal  surface  at  the  bottom  of  the  cell,  the  Nusselt 
number  is  also  given  by  Eq.  (13)  [11]: 

Nu  =  0.54  Ral/A.  (13) 


For  vertical  surfaces,  the  corresponding  Nusselt  number  is 
given  by  (14),  considering  a  laminar  flow  [12]: 

Nu  =  0.59  Ral/4.  (14) 


The  coefficient  of  heat  transfer  by  convection  permitting  to 
estimate  the  conductances  is  obtained  from  Eq.  (15): 


Nu  x 
h  =  - 

Afl 


(15) 


with  v  the  characteristic  length  of  the  considered  surface. 


temperature  of  the  cooling  system  (°C) 


Fig.  4.  Temperatures  as  functions  of  water  temperature.  Top  surface,  without 
load,  3  cells. 


cell  with  the  load  to  freeze  all  the  other  parameters  and  especially 
the  high  influence  of  the  water  temperature.  At  the  opposite, 
with  no  load,  the  modification  of  the  water  temperature  consid¬ 
ered  as  the  functioning  temperature  of  the  cell  permits  to  define 
the  unknown  internal  thermophysical  parameters  as  the  thermal 
conductivities  or  the  thermal  contact  resistances  necessary  for 
the  modeling.  For  these  six  modeling,  we  describe  the  external 
temperatures  so  that  the  3D  internal  and  external  heat  fluxes. 
Two  PEMFC  of  3  and  20  cells  with  respective  power  of  150 
and  500  W  are  studied.  Fig.  4  gives  the  results  of  measured  and 
computed  temperatures  of  tests  without  load  for  different  input 
temperatures  of  the  water.  Curve  in  configuration  1  gives  the 
measured  temperatures  obtained  thanks  to  the  infrared  camera 
on  the  isothermal  surfaces  of  the  stack  (3  cells)  with  temperatures 
of  the  cooling  water  given  at  40,  50  and  60  °C.  Two  other  curves 
give  the  computed  temperatures  for  respectively  configurations 
2  and  3. 

The  three  profiles  are  rectilinear.  Differences  of  measured 
and  computed  temperatures  on  the  top  surface  of  the  stack  do 
not  exceed  2.5  °C  for  the  temperature  of  the  water  at  50  °C. 

Fig.  5  gives  the  same  computed  and  measured  temperatures 
with  loads  at  0.3,  0.4  and  0.5  A  cm-2.  The  differences  between 
the  measured  and  the  computed  values  are  then  smaller  and  do 
not  exceed  1.5  °C.  Nu  obtained  from  Eq.  (13)  are  28. 

Fig.  6  gives  the  same  results  as  Fig.  5  for  the  20  cells  stack. 
Differences  of  the  measured  and  the  computed  values  are  similar 
as  for  the  3  cells  stack  (about  1.3  °C  for  a  temperature  of  the 
water  of  50  °C.  The  temperature  at  the  top  surface  decreases  as 


6.  Results 

Solved  equations  are  given  thanks  to  expression  (2).  Expres¬ 
sions  of  Eqs.  (2a)-(2c)  are  obtained  with  precedent  dimension¬ 
less  numbers. 

Three  loaded  experimental  tests  (0.3, 0.4  and  0.5  A  cm-2)  are 
presented  with  their  corresponding  modeling  for  the  output  tem¬ 
perature  regulated  at  50  °C.  Other  experimental  phases  with  no 
load  are  realized  for  40,  50  and  60  °C  as  output  water  tempera¬ 
tures.  The  reason  of  the  constant  temperature  of  the  water  for  the 
tests  with  load  is  to  particularize  the  internal  power  sources  of  the 
cell  as  function  of  the  load.  Indeed,  it  is  fundamental  for  under¬ 
standing  the  electrical,  chemical  and  thermal  behaviors  of  the 


current  density  (A/cm2) 


Fig.  5.  Temperature  on  the  top  side,  with  load,  3  cells. 
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Table  3 

Characteristics  of  the  flow  in  the  main  duct 


40  °C 

50  °C 

60  °C 

Without  load 

3  cells 

AT  (°C) 

1 

0.6 

0.9 

Llow  rate  (lmin-1) 

1.06 

1.09 

0.91 

0.3  A  cm-2 

0.4  A  cm-2 

0.5  A  cm-2 

With  load  (T=50°C) 

3  cells 

AT  (°C) 

0.2 

-0.3 

-0.6 

Llow  rate  (lmin-1) 

1.02 

1.00 

1.08 

20  cells 

AT  (°C) 

0.9 

-0.1 

-1.3 

Llow  rate  (lmin-1) 

2.41 

2.47 

2.45 

the  current  density  increases.  Indeed,  the  operating  temperature 
is  controlled  with  the  output  temperature  of  the  cooling  system. 
The  heat  released  by  the  electrochemical  reaction  increases  with 
the  current  density.  The  input  temperature  of  the  cooling  system 
decreases  as  considering  constant  the  output  temperature.  As  the 
current  density  increases,  the  electrochemical  reaction  produces 
more  heat.  Temperatures  of  the  stack  tend  to  become  higher 
as  the  temperature  of  the  feeding  system.  This  effect  is  more 
important  for  the  20  cells  stack  where  the  heat  released  by  the 
electrochemical  reaction  is  more  important  than  for  the  3  cells 
stack. 

Finally,  after  having  integrated  modeled  thermal  fluxes  of  the 
cell  on  the  entire  PEMFC,  the  balance  of  power  of  the  consid¬ 
ered  cell  is  done.  To  reach  that  result,  it  is  necessary  first  of  all  to 
estimate  experimentally  the  heat  fluxes  transferred  through  the 
water  duct.  The  heat  fluxes  transferred  to  the  other  fluids  are  con¬ 
sidered  as  negligible.  Indeed,  the  experimental  values  obtained 
at  the  higher  load  at  50  °C  are  respectively  0.7  and  1.8  W  for 
one  cell  (case  of  the  3  cells  stack)  for  the  hydrogen  and  the  air 
ducts.  The  sum  of  these  ones  does  not  exceed  1%  of  the  heat 
power  transferred  with  all  the  fluids.  Table  3  gives  the  experi¬ 
mental  results  with  A  Fas  the  difference  of  the  input  and  output 
temperatures  of  water. 

Table  4  is  the  power  balance  according  to  the  first  law  of 
thermodynamic  for  the  three  configurations.  The  positive  alge¬ 
braic  values  are  relative  to  gaining  of  power  in  comparison  to 


Table  4 


Energy  balance  of  the  stack 


Load  (T=50°C) 

0.3  A  cm  2 

0.4  A  cm  2 

0.5  A  cm  2 

3  cells 

Chemical  (W) 

112.96 

150.50 

187.82 

Electrical  (W) 

-60.23 

-72.89 

-85.74 

Thermal  (W) 

-52.73 

-77.61 

-102.08 

20  cells 

Chemical  (W) 

749.35 

998.64 

1253.47 

Electrical  (W) 

-401.80 

-514.24 

-640.85 

Thermal  (W) 

-347.54 

-484.40 

-612.62 

the  negative  values  for  power  losses  by  the  stack.  Table  5  gives 
the  different  computed  leak  conductances  to  balance  the  thermal 
power  sources  in  the  stacks  for  the  three  configurations. 

Table  6  gives  the  same  computed  parameter  for  the  two  stacks 
with  load,  the  stack  is  regulated  at  50  °C. 

It  can  be  seen  in  Table  6  that  the  leak  conductance  becomes 
smaller  as  the  current  density  increases.  The  relative  part  of  heat 
transfer  extracted  from  the  stack  by  the  cooling  system  increases. 
The  input  temperature  of  the  cooling  system  increases  as  the 
current  density  is  smaller,  being  also  more  efficient.  The  leak 
conductances  for  the  configuration  3  are  then  smaller  compared 
to  the  other  configurations.  The  taking  into  account  of  the  bound¬ 
ary  conditions  defined  in  configuration  3  between  two  adjoining 
cells,  creates  a  heat  flux  on  the  axes  O ,  z  which  reduces  the  part 
of  the  heat  lost  with  the  other  kinds  of  heat  transfer. 

Fig.  7  quantifies  the  heat  fluxes  of  one  computed  cell  from 
the  cathode  bipolar  plate  to  the  anode  bipolar  plate  of  the  next 
(Active)  cell.  These  thermal  fluxes  are  expressed  as  functions  of 
the  current  densities  produced  by  the  cell. 

The  two  curves  are  relative  to  the  two  stacks.  Configura¬ 
tions  1  and  2  are  finally  similar  considering  that  they  do  not 
allow  heat  losses  by  conduction  in  direction  O ,  z  (homogeneous 
Neumann  conditions  and  cyclical  conditions).  For  that  reason, 


Table  5 

Leak  conductances,  without  load,  3  cells 


40  °C 

50°C 

60  °C 

Configuration  1  (WK-1) 

2.18 

0.86 

0.71 

Configuration  2  (W  K- 1 ) 

2.19 

0.87 

0.73 

Configuration  3  (W  K- 1 ) 

1.75 

0.44 

0.31 

Table  6 

Leak  conductances,  with  load,  3  and  20  cells 


0.3  A  cm  2 

0.4  A  cm  2 

0.5  A  cm  2 

3  cells 

Configuration  1  (WK-1) 

0.88 

0.80 

0.78 

Configuration  2  ( W  K- 1 ) 

0.89 

0.80 

0.78 

Configuration  3  ( W  K- 1 ) 

0.44 

0.35 

0.32 

20  cells 

Configuration  1(WK-1) 

7.09 

6.79 

5.95 

Configuration  2  ( W  K- 1 ) 

7.09 

6.79 

5.95 

Configuration  3  ( W  K- 1 ) 

7.02 

6.72 

5.88 

Lig.  6.  Temperature  on  the  top  side,  loaded  tests,  20  cells. 
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current  density  (A/cm2) 

Fig.  7.  Heat  flux  from  the  cathode  bipolar  plate  to  the  anode  bipolar  plate,  in 
configurations  1  and  2. 

configuration  3  is  the  more  realistic  case,  but  it  is  very  difficult 
to  quantify  with  a  very  good  accuracy. 

Figs.  8  and  9  show  respectively  the  heat  transfer  in  the  direc¬ 
tion  O ,  z  for  one  computed  cell  for  both  stacks.  The  heat  fluxes 
are  lost  for  the  system  and  dispelled  to  the  external  medium. 
They  are  expressed  as  a  function  of  the  current  densities  (con¬ 
figuration  3).  The  difference  of  the  influence  of  the  cooling  fluid 
can  be  seen  on  both  figures  through  the  values  of  the  flux.  It  is 
shown  in  that  way  the  high  influence  of  the  cooling  system  for 
the  stack  of  20  cells  (as  shown  in  following  Table  8).  Thermal 
conductances  between  the  stack  and  the  external  medium  are 
defined  as  a  function  of  the  number  of  cells  for  each  stack.  Con¬ 
sequently,  the  axial  heat  flux  is  higher  for  the  3  cells  stack  than 
for  the  20  cells  stack. 


Fig.  8.  Heat  flux  to  external  medium  in  “z”  axis,  case  3,  with  load,  3  cells. 


Table  7 


Heat  energy  extracted  by  the  cooling  system 


0.3  A  cm  2 

0.4  A  cm  2 

0.5  A  cm  2 

3  cells  (%) 

22.2 

-24.6 

-43.4 

20  cells  (%) 

42.8 

-2.7 

-35.4 

Table  8 

Heat  energy  extracted  by  the  cooling  system 

0.3  A  cm  2 

0.4  A  cm  2 

0.5  A  cm  2 

3  cells 

Configuration  1 

-86.3 

-88.7 

-88.6 

Configuration  2 

-88.9 

-90.1 

-89.5 

Configuration  3 

85.7 

29.7 

2.5 

20  cells 

Configuration  1 

-88.6 

-90.2 

-90.2 

Configuration  2 

-90.4 

-91.6 

-90.9 

Configuration  3 

-62.3 

-71.2 

-75.2 

Tables  7  and  8  compare  and  quantify  the  relative  influence  of 
the  cooling  system  at  different  current  densities  for  one  cell  and 
all  the  cells  of  the  two  stacks.  As  shown  in  Table  7,  both  stacks 
are  not  able  to  reach  the  operating  temperature  at  0.3  A  cm-2. 
The  “cooling”  system  is  then  a  heater  system.  Table  8  describes 
the  dissipation  way  of  the  heat.  Configurations  1  and  2  in  both 
stacks  are  similar. 

When  heat  cannot  be  dispelled  in  the  axial  axis  O ,  z  because 
the  modeling  does  not  take  into  account  this  possibility  (config¬ 
urations  1  and  2),  the  cooling  channels  are  the  main  processes 
to  extract  the  heat.  Heat  transfers  by  conduction  along  axes  O ,  v 
and  (9,  y  are  not  very  important,  particularly  for  the  configura¬ 
tion  2  with  the  computed  wall  temperatures.  The  configuration  3 
takes  into  account  the  heat  flux  on  axis  (9,  z.  For  the  3  cells  stack, 
the  cooling  channels  are  in  fact  some  heaters  for  all  values  of 
the  current  density.  It  is  not  realistic.  For  the  20  cells  stack,  even 
in  the  configuration  3,  the  heat  flux  is  mainly  removed  by  the 
cooling  channels.  The  thermal  conductance  in  (9,  z  axis  is  not  a 
linear  function  of  the  number  of  cells  and  the  effective  thermal 
conductance  is  smaller  than  the  physical  thermal  conductance 
computed  from  the  different  components  of  one  cell. 

In  both  Tables  7  and  8,  the  positive  values  correspond  to 
contribution  of  power  and  negative  values  to  the  heat  losses 
when  considering  the  stack  as  a  thermodynamical  system. 


Fig.  9.  Heat  flux  to  external  medium  in  “z”  axis,  configuration  3,  with  load,  20 
cells. 


7.  Discussion 

Heat  transfers  in  the  stack  change  with  the  operating  con¬ 
ditions.  When  the  stack  operates  without  load  (3  cells  stack), 
the  feeding  system  is  used  to  reach  the  operating  temperature. 
For  a  water  output  temperature  regulated  at  40  °C,  the  heat  flux 
between  two  cells  is  negligible.  Indeed,  the  difference  of  tem¬ 
peratures  between  the  wall  and  the  external  air  is  about  15  °C. 
The  different  cells  also  operate  under  adiabatic  conditions  (con¬ 
figuration  2).  As  the  temperature  of  the  water  increases  the  heat 
fluxes  in  the  axis  O ,  z  follow  the  same  evolution.  The  stack 
operates  also  according  to  the  third  configuration.  For  a  cur¬ 
rent  density  of  0.3  A  cm-2,  the  heat  quantity  created  by  the 
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electrochemical  reaction  stays  insufficient  to  maintain  constant 
the  operating  temperature.  The  heat  power  released  by  the  elec¬ 
trochemical  reaction  is  also  extracted  by  conduction.  The  feed 
system  then  behaves  as  a  heater  system.  The  boundary  conditions 
of  the  cells  stay  near  to  those  modeled  in  configuration  3.  With 
load,  the  heat  released  by  the  electrochemical  reaction  increases 
with  the  current  density.  For  a  current  density  of  0.4  A  cm-2 
the  feed  system  extracts  the  heat  generation  for  the  three  cells 
stack  whereas  the  equilibrium  of  fluxes  is  quite  reached  for  the 
twenty  cells  stack  (Table  7).  As  the  current  density  goes  up,  it 
can  be  seen  a  difference  of  temperatures  of  the  system  between 
the  configurations  3  and  2.  It  means  that  the  system  moves  to 
an  adiabatic  mode.  For  the  20  cells  stack,  it  can  be  observed 
a  small  difference  between  configurations  2  and  3.  Indeed,  the 
thermal  conductivity  is  supposed  inversely  proportional  to  the 
number  of  the  cells.  It  is  estimated  to  about  0.02  W  m-1  K-1  for 
the  20  cells  stack  whereas  it  was  about  0. 15  W  m  1  K  1  for  the 
3  cells  stack.  Nevertheless,  the  wall  temperatures  computed  are 
lower  than  the  measured  temperatures.  The  internal  heat  fluxes 
distributions  are  different  between  the  two  stacks  because  of  the 
important  differences  of  the  geometries  at  the  origin  of  different 
internal  thermo  physical  parameters  (convective  coefficients  so 
that  equivalent  thermal  conductivities). 

8.  Conclusion 

To  realize  the  thermal  modeling  of  a  PEMFC,  it  is  necessary 
to  predict  internal  temperatures  reached  for  different  function¬ 
ing  configurations.  Finally,  the  global  efficiency  of  the  cell  with 
its  feed  system  can  be  estimated.  To  validate  the  model  devel¬ 
oped  for  one  cell  (instead  of  3  or  20),  it  has  been  necessary  to 
obtain  the  convergence  of  computed  and  simulated  temperatures 
at  the  output  of  the  cooling  water.  The  approach  is  particularly 
complicated  in  reason  of  the  relative  absence  of  experimental 
measurements  on  the  stack.  The  validation  of  the  modeling  tool 
is  realized  for  different  input  temperatures  of  water  (without 
load)  and  with  different  loads  ranging  from  0.3  to  0.5  A  cm-2. 
The  output  temperature  of  the  cooling  channel  is  regulated  and 
also  considered  as  constant  at  50  °C.  The  validated  tool  has  also 


permit  to  estimate  the  internal  and  external  ways  of  heat  transfer. 
Studies  of  one  cell  and  for  the  entire  stack  (3  and  20  cells)  and 
their  integration  in  a  global  thermal  behavior  model  will  then 
permit  to  estimate  the  final  efficiency  of  the  stack  with  the  feed 
system.  Next  step,  according  to  the  important  losses  and  non 
utilized  external  heat  fluxes,  will  consist  in  the  development  of 
the  utilization  of  this  power  on  the  test  bench  to  optimize  the 
efficiency  of  the  stack  with  its  feed  system.  For  this  it  is  neces¬ 
sary  to  modelize  more  than  one  cell  in  order  to  take  into  account 
the  influence  of  the  neighboring  cells  to  consider  the  heat  flux 
distribution. 
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